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Summary 


In this paper we first introduce the least-squares (L 2 ) finite element method for two- 
dimensional steady-state pure convection problems with smooth solutions. We prove that 
the L 2 method has the same stability estimate as the original equation, that is, the L 2 
method has better control of the streamline derivative. Numerical convergence rates are 
given to show that the L 2 method is almost optimal. Then we use this L 2 method as a 
framework to develop an iteratively reweighted L 2 finite element method to obtain a least 
absolute residual [L x ) solution for problems with discontinuous solutions. This L x finite 
element method produces a non-oscillatory, non-diffusive and highly accurate numerical 
solution that has a sharp discontinuity in one element on both coarse and fine meshes. We 
also devise a robust reweigting strategy to obtain the L x solution in a few iterations. A 
number of examples solved by using triangle and bilinear elements are presented. 
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1. Introduction 


In this paper we introduce and test numerically the Li finite element method for the 
solution of two-dimensional steady-state pure convection problems. This new method is 
designed to obtain accurate non-oscillatory discontinuous solutions. We shall consider the 
following steady-state boundary value problem: 

Uj) — 0 in fi, (l-°) 

u = g on T_, (1-^) 

where fi is a bounded convex domain in 9J 2 with boundary r, u — u{x,y ) is the dependent 
variable (e.g., the concentration), f3 = (& ,/? 2 ) is a constant vector with \f}\ = 1, Up = fi -Vu 
denotes the derivative in (3 direction, and g is the given data on the inflow boundary T_ 
defined by 

r_ = {(*, y) € r : n(i, y) • (3 < 0}, 

in which n is the outward unit normal to T at point (x, y) E T. The problem (1) is purely 
hyperbolic. The characteristics of the problem (1) are the straight lines parallel to /?. 
The analytic solution of problem (l) is very simple. The solution is a constant along a 
characteristic. The value of this constant is equal to the given value of g at the intersection 
of this characteristic and the inflow boundary. However, the solution is discontinuous with 
a jump across a characteristic, if the boundary data g is discontinuous. This creates the 
greatest challenge to numerical solutions. 

Commonly used numerical methods for hyperbolic problems are of the following 
types(see e.g., Fletcher[lO], Hirsch[I3], Johnson[20] and Pironneau[29]): method of charac- 
teristics, finite difference and finite element methods. In principle the method of character- 
istics is very good, but it is rather cumbersome in practice. Usually one uses finite difference 
and finite element methods based on a mesh, which is not adapted to fit the characteristics 
of the particular problem. In such a case, if the exact solution has a jump discontinuity 
across a characteristic, all conventional finite difference and finite element methods will 
produce approximate solutions which either oscillate or smear out a sharp front. Finding 
accurate approximations of the discontinuous solutions of hyperbolic equations has been a 
persistent difficult task in modern numerical mathematics and computational physics. 

One research direction towards better resolution around discontinuities is to use an 
adaptive h-refinement strategy, such as that extensively investigated by Oden and his 
colleagues [2 8]. However, the data structure and the programming of h-refinement are 
complicated, especially for three-dimensional problems. 

Impressive sharp discontinuities may be obtained by using filter methods[5,22,9]. For 
example, one may use TVD-type finite difference schemes[ll,12] to get a good approxima- 
tion, then use filter methods to improve the resolution. 
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Another potential way is to use conventional finite difference or finite element meth- 
ods to get approximate solutions, then apply imaging processing techniques to detect the 
locations of discontinuities [30,31]. 

The Li procedure for non-oscillatory solutions, first proposed by Lavery in [23,24], 
is a rather different approach. The L x idea can be explained as follows. In the usual L 2 
curve fitting, the L 2 procedure does its best in a sense of least-squares of the residual to 
make the curve pass through or by all of the data. If the data are smooth, the L 2 fitting 
leads to a very good approximation. However, if the data contain abrupt changes, the L 2 
procedure will produce an oscillatory and diffusive curve around sharp changes. In such a 
case, the trouble comes from the fact that the L 2 fitting makes the use of individual datum 
equally important. The tendency of L x fitting is to give up the outliers in the data and to 
require the remaining data be satisfied exactly. Therefore, the L x fitting is the choice for 
discontinuous functions. The same thing happens in the L 2 and L x solutions of discretized 
hyperbolic equations. The L x capacity translates into a capacity to permit the equation 
in the “shocked” cell (in which the discretized scheme does not hold) not to be satisfied 
while requiring that the remaining equations be satisfied exactly. The L x solutions are 
non-oscillatory, highly accurate and right up to the edge of the discontinuity. 

The key point of L x procedure is how to get overdetermined discretized systems. The 
standard finite difference and finite volume methods lead to determined linear systems. In 
order to get an overdetermined linear algebraic systems, one must rely on non-traditional 
tricks, such as, artificially adding the extra boundary conditions on the outflow boundary 
for two-dimensional pure convection problems [25], or gradually adding a smaller viscous 
term for one-dimensional Burgers’ equation[23]. Another big difficulty associated with the 
usual L x procedure is that the linear programming algorithm of Barrodale and Roberts] 1] 
is very expensive. It requires at least 0(n 4 ) and perhaps as many as 0{n 6 ) operations, 
here n is the number of grids in each axis. This excludes the possibility of practical use of 
the usual L x method. 

In this paper we first introduce the least-squares (L 2 ) finite element method for two- 
dimensional steady-state pure convection problems with smooth solutions. We prove that 
the L 2 method has the same stability estimate as the original equation, that is, the L 2 
method has better control of the streamline derivative. Numerical convergence rates are 
given to show that the L 2 method is almost optimal. Then we use this L 2 method as a 
framework to develop an L x finite element method for the solution of hyperbolic equations. 
The L 2 finite element method with numerical quadrature is equivalent to a weighted col- 
location least-squares method[4], in which at first the residual equations are collocated at 
the interior points in each element, then the algebraic system is approximately solved by 
the weighted least-squares method. The Gaussian points for calculating the element ma- 
trices in the L 2 finite element method correspond to the collocation points in collocation 
methods. If the order of Gaussian quadrature (or the number of quadrature points) is 
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appropriately chosen, the L 2 finite element method amounts to solving an overdetermined 
system. 

Since the L 2 finite element method produces a very good approximation to the exact 
solution, we may use the information provided by the L 2 solution to find “shocked” ele- 
ments. Then we use the L 2 method again, but this time we put a small weight for “shocked” 
elements, and repeat this procedure a few times until the solution is reached. 

The arrangement of this paper is as follows. The L 2 method and the convergence 
tests for smooth problems are presented in Section 2. In Section 3 we describe the L x 
procedure. The numerical results in Section 4 contains the Li solutions for the pure 
convection problems with discontinuities. Conclusions are drawn in Section 5. 


2. The L 2 Finite Element Method 

2.1 Preliminaries and Notations. As we have already shown in our previous papers (see 
[19] and the references therein), the L 2 finite element method is a universal method for 
the numerical solution of any type of partial differential equations. It does not matter 
whether the partial differential equations are elliptic, parabolic or hyperbolic. As long as 
the partial differential equation has a unique solution, the L 2 finite element method always 
gives a reasonably good approximate solution. The work done in this paper is a natural 
extension of our L 2 method. 

The problem (1) can be taken as a time-dependent problem, if we consider one space 
coordinate as a time-like coordinate. Then we may use the implicit time-marching L 2 finite 
element method introduced in [3] to get an approximate solution. The time marching is 
necessary for the Euler equations in aerodynamics [17, 18], since the Euler equations have 
nonunique solutions. The time-marching L 2 finite element method implicitly introduces 
an artificial dissipation to exclude the solutions with expansion shocks. For the linear 
hyperbolic problem (1), the time-marching is not necessary, because it has a unique so- 
lution, continuous or discontinuous. For this reason, we rather treat the problem (1) as 
two-dimensional, and directly apply the L 2 finite element method to attack it. The general 
formulation of L 2 finite element methods for first-order partial differential equations can 
be found in [19]. 

In order to compare the L 2 finite element method with other finite element methods, 
we would like to discuss more details here. Let us consider the following more general 
linear hyperbolic equation: 


Up + u — f in fl, 

(2.a.) 

u = g on T_ , 

(2.6) 
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where / is a given source function. Without loss of generality we assume that the boundary 
data g is zero. 

Throughout this paper, we use the following notations. L 2 (fl) denotes the space of 
square-integrable functions defined on fl with the inner product 

u,v € l 2 (g), 

u e L 2 (n). 

Li (fl) denotes the space of functions defined on fJ, whose absolute values are integrable 
with the norm 

HU, = / |u|dfl tt G Li(fl). 

J n 

H r (0) denotes the Sobolev space of functions with square-integrable derivatives of order 
up to r. || . || r denotes the usual norm for H T (Q). We also use the following notations: 


and the norm 


(u,v) = / 
J n 


uvdCl 


u 


= («,«) 


< 


u,w >= J 

L 


< U,W >+■ 


uwn • fids, 
uwn ■ (ids, 


ju| r =(/ « 2 |n • /3\ds)* , 


where 


r+ ={(x,y)Gr:n(a:,y)-/?>0}. 
We note that by Green’s formula 


{u 0 ,w) =< u,w > -(u,w 0 ). 


Further we also define the function space 

S = {u G H l (Q) : u = 0 on T_ }, 

and the corresponding finite element subspace S h , i.e. S h is the space of continuous 
piecewise polynomial functions of degree k. Here the parameter h represents the maximal 
diameter of the elements. By the finite element interpolation theory[7,27] we have: Given 
a function u £ H k + l (fl), there exists an interpolant u fc € S h such that 

||u-u h || <c/i* + 1 |M| fc+1 , (3.a) 
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|| Vu- V« fc || <c/i fc ||u|| fc+1 . ' (3.6) 

2.2 The Standard Galerkin Method. Now let us look at the following standard Galerkin 
method for the problem (2) (see e.g., Johnson[20|): Find u h G S h such that 

(u h 0 +u h ,w h ) = (f,w h ) Vw h G S h . (4) 

We define the error e = u —u h , then the error estimate for the standard Galerkin method 

ll e ll 4" |e|r < ch k ||u||fc+ 1 , (5) 

which is one order lower than that for elliptic and parabolic problems. Furthermore, in 
the continuous problem (2), we have the following stability estimate: 

INI + 11% II + Mr < c||/||, (6) 

But in the standard Galerkin method the stability estimate is 

||«‘|| + |«‘lr <c|l/ll. P) 


which has no control of ||u£ ||. 

2.3 The SUPG Method. In order to get better accuracy and stability, methods of upwinding 
type have appeared, see e.g., papers by Dendy[8], Wahlbin[32], Christie, Griffths, Mitchell 
and Zienkiewicz[6], Hughes and Brooks[l5], Johnson, Navert and Pitkaranta[2l], Morton 
and Parrot[26|. Below we shall look at the streamline upwinding Petrov-Galerkin method 
(SUPG)(Hughes[14|) or the streamline diffusion method(Johnson(20]): Find u h € Sh such 
that 

(lip + u h ,w h + hw h ) = (f,w h + hw h ) Vw h es h . (8) 

Johnson and his colleagues have derived the error estimate: 

+ + <«*‘ +} IMU + .. ( 9 ) 

which is near optimal. However, in SUPG the corresponding stability estimate is 

IN ll + vNNii + Nlr ^ c ii y iu (^-®) 

which means that the streamline derivative is less controlled. Another disadvantage of the 
SUPG in practical calculation is that the stiffness matrix is non-symmetric which makes 
the solution of large-scale problems very difficult. 

2.4 The L 2 Method. Now let us introduce the L 2 finite element method. We assume 
that / G L 2 (n). For an arbitrary trial function v G 5, we define the residual function 
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R = Vp + v — f. The L 2 method is based on minimizing the residual function in a least- 
squares sense. We construct the least-squares functional: 

I ( v ) = ll^ll 2 = bn + v - /|| 2 = {v 0 + v - f,v p + v - /). (11) 

The L 2 method reads: Find u E S such that 

J(u) < I(v) Vu € S. 

Taking variation of I with respect to v, and setting SI = 0 and Sv = w, lead to the L 2 
weak statement: Find u E S such that 

b(u,w) = l(w) Vw E S, (12) 

where b(u,w ) = ( u 0 + u,w 0 + w ) and l(w) — + w). The corresponding L 2 finite 

element method has the following form: Find u h E S h such that 

b[u h ,w h ) = l{w h ) Vw h E S h . (13) 

Let us now turn to the error estimate. Since we can replace w in (12) by w h , we have 

b(u,w h ) = l(w h ) Vw h E S h . (14) 

By subtracting (13) from (14) we get the following orthogonality for the error e: 

b(e,w h ) = 0. 

Let u E S h be the interpolant of u satisfying (3) and write p = u — u h and 6 = u h — u h so 
that e — p + 6. Then we have 

\\ e i) + e l| 2 = &(e, e) = 6(e, p) + b(e , 6) = b(e, p) 

< \\ e 0 + *1111 P0 +HI» 

or 

\\ e 0 + g ll ^ II Pfi + P|| < \\Pp || + ||/>||- 
Recalling (3) we obtain the error estimate: 

11^/? + e|| < ch fc ||u|| fc+i . (15) 

Since the residual of approximate solution R h = u h p + u h — / = e 0 + e, (15) is also the 
residual estimate: 

ii^ii<^w +i , (i6) 
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which means that the residual estimate is optimal. Using Green’s formula and the bound- 
ary condition, (16) can be rewritten as 

(l|e|| 3 + l|e»ir+<e. e > + )i<cA*||„|U +1 , (17) 

which shows that the error estimate for e 0 is optimal, but the error estimate for e is one 
order lower than optimal. Although in numerical tests (see below) we have observed that 
the accuracy of the L 2 method is higher than the fcth order, it is still an open question for 
getting a better theoretical error estimate in general. 

By taking w h = u h in (13), we can obtain the stability estimate: 

IKII + Kll + I^lr <c||/||. (18) 

This estimate is the same as the above estimate (6) for the continuous problem. It means 
that the L 2 method has better control of the streamline derivative. We also note that the 
bilinear form in (13) is symmetric, therefore the matrix of the resulting algebraic system 
is symmetric and positive definite. This is a very important advantage of the L 2 method 
over other methods in practice. 

2.5 Numerical Experiments of the L 2 Method. We chose the following model problem: 

du du . . . _ , 

- — b — = «mi + y) m U, (19. a) 

ox oy 

u — 0 on T_ , (19.6) 

where Q — {(z,y) 6 SR 2 : 0 < x < 1, 0 < y < 1} is the unit square, and T_ = {(x, y) 6 T : 
x = 0 or y = 0}, in which T is the boundary. This problem has a smooth exact solution 
u — sin(x)sin(y). 

We have tested the bilinear element with uniform meshes. At first the one-point 
Gaussian quadrature was used for calculating the stiffness matrices. In the case of one- 
point quadrature, the L 2 method is equivalent to the collocation least-squares method 
with one collocation point at the center of each element. It is easy to check that in such a 
collocation method, the number of discretized algebraic equations “ nequ ” is equal to the 
number of unknowns “nelem 2 ”, here “ nelem ” is the number of elements. In other words, 
in such a case we solve a determined system. Therefore, there is no difference between the 
L 2 solution and the direct collocation solution. The numerical result for convergence rate 
is shown in Figure 1. The optimal rate, i.e. J|e|| < ch 2 , is observed. 

Also we would like to mention that the L 2 method with one-point quadrature is 
equivalent to the central finite difference scheme. Therefore, the ^ 2 -optimality may be 
derived by using the finite difference theory. 
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The numerical rate of convergence with the 2x2 Gauss rule is also included in Figure 
1. In this case, the L 2 method solves an overdetermined system. The convergence rate is 
around 0(h 1 75 ), which is near optimal. Here, more theoretical study is needed. 

We also did the numerical tests with specified extra boundary conditions on the out- 
flow boundary T + = {(x, y) € T : x = 1 or y — l}. In this case, the L 2 method with the 
2x2 Gauss rule gives the optimal rate of convergence ||e|| < ch 2 (Figure 1). 


3. The Finite Element Method 

If the solution is non-smooth, the above L 2 method still performs quite well. Of course, 
the argument about the error estimates does not hold. As expected, the L 2 method smears 
out the jump across a characteristic. As we pointed out in Section 1, the trouble comes from 
“shocked” elements, where the direction derivative across the jump approaches infinity and 
the discretized equation is not valid. But the usual L 2 method does not recognize them, 
and just equally treats “shocked” and “smooth” elements. This is the reason that we 
would like to use the Li idea to suppress the interference of “shocked” elements. 

We still consider the problem (2). We want to minimize the L x norm of the residual: 

\\R\\ Ll = / \u 0 +u-f\dn. ( 20 ) 

Jo 

Obviously, we should work on the discretized version of (20) by using the finite element 
method: i.e. Find the minimizer u h of 

nelem n gaut 

ll«*IU. = E E «*IW(6.'J|)I, ( 21 ) 

J = 1 (= 1 

in which 

Ri = + u' l (6,»?i) - f{Zt,Vi), ( 22 ) 

where R t stands for the residual at each Gaussian point, ngaus denotes the number of 
Gaussian points, iuj is the Gaussian weighting, |J| is the determinant of the Jacobian 
matrix, and (Ci,Vi) is the local coordinates of Gaussian points. As usual, u h can be 
expressed by the shape functions and the nodal values: 

n n ode 

u h {t,ri) = £ 

m — 1 

where “ nnode ” is the number of nodes in an element, are the shape functions, and 
U m are the nodal values. In order to make the problem (21) meaningful, we must have an 
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overdetermined algebraic system. It can be realized simply by appropriately choosing the 
number of Gaussian points. 

Now the task becomes finding the solution which has the least absolute residual. It is 
well known that any L x problem can be transferred into a linear programming probIem[2]. 
Then we may use a linear programming algorithm to find the L 1 solution. But right now 
this type of algorithm is too time-consuming. We are still waiting for a fast algorithm 
which can at least take the advantage of sparse finite element matrix. 

Fortunately, we can use another approach, e.g., an iteratively reweigh te least-squares 
method(IRZr 2 )[2], which is based on repeatedly solving a weighted least-squares problem: 
Find the minimizer u h of 




nelem ngaus 

£ E 7i)i, 

y=i /=i 


(23) 


in which 

(24) 

where W t denotes the weight set, which in turn depends on the information of the previous 
step. The IRX 2 would begin with the initial weight set W t = 1. This first step is nothing but 
the L 2 method introduced in Section 2. The result of this L 2 method determines a new set 
of weights by (24). In the second iteration, the residual |/ij| is larger in “shocked” elements. 
Thus the weight W, for “shocked” elements is smaller, and their inference becomes less 
important. This procedure is repeated until ||u!* — u h II is small. 

1 II CUfrviil pf 6 V IOU I H 

Our numerical experiments reveal that the above method converges extremely slowly. 
The problem is that the difference between the residuals of “shocked” and their neighboring 
elements in the first L 2 solution is not significant enough. This difficulty can be overcome 
simply by using the following weights: 


W t 


R, 


previou 8 


w, = 


\R, 


previous 


(25) 


It corresponds to additionally increasing the importance of “smooth” elements and reduc- 
ing the inference of “shocked” elements. This is reasonable, since the theory of L, fitting(2] 
tells us that the L x procedure eliminates completely the equations, which will have nonzero 
residuals, from the system. This trick is usable, also because our non-weighted L 2 method 
is good enough to locate the “shocked” elements. That is, in the results of our L 2 method, 
the absolute value of the residuals in “shocked” elements is always greater than that in 
other elements. 
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We may further simplify the procedure by using another simple and reliable “shock” 
indicator-the variation of nodal values in each element- instead of the residual. The 
variation is defined as 


nnode 

V =Y1 U 0 =U nnode . (26) 

m = 1 

The advantage of using the variation as a “shock” indicator is as follows: Once the jump 
in the boundary data is given, we may know the exact values of the variation in “shocked” 
elements in advance. There are only a few possible values, which depend only on the type 
of finite element and are independent of the shape and size of the particular element, and 
have no relation with the location of quadrature points. 

The implementation of this L\ method is really straightforward. If an L 2 finite element 
code is already available, it needs only a few additional lines of FORTRAN statements. 


4. Numerical Results of the L 1 Method 
We consider the following problem: 

^ + ian(35*)^ = 0 in U, (27. a) 

ox ay 

where 0 = {(x,y) € 9J 2 : 0 < x < 1, 0 < y < 1} is the unit square with the boundary T. 
The inflow boundary conditions are 

u — 2 on Tj = {(x,y) £ T : x = 0}, (27 .6) 

u = 1 on r 2 = {(x,y) £ T : x > h and y = 0}, (27. c) 

in which h is a positive constant less than 1 ( h will be a mesh length). Equations (27) 
represent uniform flow along straight lines inclined at an angle of 35° with respect to the 
x-axis. In this case, any straight line, which is between and parallel to y = xtan( 35°) and 
y = (x — h) tan (35°), could be considered as the location of discontinuity. For example, we 
may write the solution of (27) as 

u = 2 on and above the line y — (x — |)tan( 35°), 

u = 1 below the line y = (x — |)tan(35°). 

The jump discontinuity occurs along the line y = (x — |)fan.(35°). 
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The boundary conditions (27.6) and (27. c) can be transferred into the source term in 
the equation (27. a). Therefore, the formulation of the L 2 method described in Section 2 
can be directly applied to the problem (27). 

The computational results presented in this paper were obtained in double precision 
on our PC-386. A direct solver with variable band-width was used to obtain the solution 
of linear algebraic equations. The computing time will be significantly cut by using the 
preconditioned conjugate gradient method[l6], since the L 2 solution is already close to L x 
solution and the final iteration is often just for correcting one or two nodal values which 
have not yet reached 15 digit accuracy. 

4.1 Linear Triangle Element. At first numerical experiments were carried out for the 
problem (27) using linear triangle elements on uniform meshes with n — 5, 15. Here n is the 
number of grids in each coordinate. For triangle elements, we use the one-point Gaussian 
quadrature. Since there are 2n 2 elements, it corresponds to having 2n 2 equations. Since 
there are (n+ l) 2 nodal values and (2n-f 1) boundary conditions, the number of unknowns 
is (n + l) 2 — (2n + 1) = n 2 . That is, the number of equations is double the number of 
unknowns. Therefore, the L 2 method amounts to solving an overdetermined system. It 
does not make sense to take more quadrature points, because in a linear triangle element 
du h j dx and du h jdy are constants, and thus the residuals at different points are the same. 

The L 2 results for n = 5 (50 triangle elements) are listed in Table 1. The numbers 
in Table 1 are the nodal values. Because the mesh is very coarse, the jump discontinuity 
is smeared severely. Starting from this bad L 2 solution, after 4 iterations of IRL 2 , we 
obtained the perfect L x solution listed in Table 2. This solution has absolutely accurate 15 
digits. Here we should note that a double precision (8 bytes or 64 bits) real number in a 
computer can only represent a decimal number with 15 digits. This solution has completely 
no oscillation and no diffusion. The transition over the discontinuity is accurately located 
in the vicinity of the line y = (x — j)tan(35°), which has a width of hsin( 35°) in a sense 
discussed above, and is accomplished in just one element. 

The L 2 solution for n = 15 (450 triangle elements) is given in Table 3. This solution 
is diffused, and slightly oscillatory around the jump. Starting from this L 2 solution, the 
accurate solution is obtained after 4 iterations of IRL 2 (see Table 4). The L x solution 
again has 15 digit accuracy. Because of the limitation of page size, we give the nodal values 
with only 8 digits in Table 4. 

4.2 Bilinear Element. Numerical experiments were also carried out for the problem (27) 
using bilinear elements on uniform meshes with n = 5,15,40, and 80. For bilinear elements 
we use the 2x2 quadrature. In each element we may write the finite element approximation 
of u as a bilinear function: 


u h (x, y) = a + bx + cy + dxy. 
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Thus the residual is 


R h = + ian(35°) — = b + tan{35°)c + d{y + tan(ZB°)x) t 

ox ay 

which means that the four discretized equations at four Gaussian points are independent 
in this case. (If the flow inclines at an angle of 45° or 135° with respect to the x-axis, we 
have three independent equations, because the location of Gaussian points is symmetric.) 
All together we have 4n 2 equations and n 2 unknown nodal values. Therefore, we deal with 
an overdetermined system. 

The L 2 solution of a coarse mesh with 5x5 bilinear elements is listed in Table 5. 
Starting from this rough solution, after 4 iterations we obtained the L x solution listed in 
Table 6. We again observed a crisp computational jump in one element. This solution is 
perfectly non-oscillatory and non-diffusive. We also list the summation of four absolute 
residuals and variations in each element in Table 7 and Table 8. The large numbers (2 for 
residual and 8 for variation) indicate the “shocked” elements. 

The L 2 solution of a mesh with 15 x 15 bilinear elements is presented in Table 9, and 
the corresponding contours are given in Figure 2. This approximate solution is reasonably 
good, although the discontinuity is smeared out, and slight oscillations occur. From this 
table and this figure, we can hardly tell where the jump is located. However, after 5 
iterations, a clean L x solution is reached (Table 10 and Figure 3). This solution again has 
the correct 15 digits. The element residuals and variations are given in Table 11 and Table 
12 to contrast “shocked” elements with “smooth” elements. 

The L 2 solution of a mesh with 40 x 40 bilinear elements is illustrated in Figure 4. 
Taking this L 2 solution as an initial solution, after 8 steps of processing, we obtained the 
L 1 solution illustrated in Figure 5. No other currently available methods can produce such 
a sharp discontinuity as this. 

We also did numerical tests for meshes with up to 80 x 80 elements, combined with 
various inflow angles and different boundary conditions. All of our results are perfectly 
accurate. Because of the page limitation, we do not present these results here. 


5. Conclusions 

A new Li procedure based on the iteration of L 2 finite element method for the solution 
of pure convection problems is developed. The overdetermined algebraic system is inher- 
ently obtained by choosing an appropriate number of Gaussian points in the formation 
of element matrices. The time-consuming linear programming for solving overdetermined 
systems turns out to be not necessary. 
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This Li finite element method captures two-dimensional discontinuity in bands of 
elements that are only one element wide on both coarse and fine meshes. The solution of 
this method has no smearing and no oscillation, and has superior accuracy. The method 
is simple and robust, and can be easily extended to three-dimensional pure convection 
problems. 

We believe that the methodology developed in this paper can be transferred into many 
other areas which deal with sharp fronts such as oil reservoir simulation, weather forecast, 
and image enhancement. We have already extended this method to two-dimensional com- 
pressible flows with shocks. 
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Figure 1. Computed convergence rate for the pure convection problem. 
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Figure 4. Contours of the L 2 solution for the pure convection problem 
(40 x 40 bilinear elements) 



Figure 5. Contours of the solution for the pure convection problem 
(40 x 40 bilinear elements) 
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0. 18702E-14 

0. 17055E- 15 

0.49194E-14 

0.48678E* 14 

0.22330E- 14 

0.92255E-15 

0.49022E-14 

0.35653E- 14 

0.11563E-15 

0. 18920E- 14 

0. 18702E- 14 

0.18828E-14 

0.59523E- 16 

0.20000E+01 

0.20000E+01 

0.85381E- 15 

0.20000E+01 

0.20000E+01 

0.20000E+01 

0. 1405QE- 14 

0.20000E+01 

0.20000E+01 

0.49960E-15 

0.64103E-16 

0 . 49960E - 1 5 

Table. 7 Element Residuals of 

Li Solution for n = 

5 ( 25 Bilinear 

Elements ) 


0.71054E-14 

0.10658E-13 

0. 19540E- 13 

0. 19540E-13 

0.88818E-14 

0, 10658E- 13 

0. 19540E-13 

0.15987E-13 

0.53291E-14 

0.88818E- 14 

0.71054E- 14 

0.88818E- 14 

0.35527E- 14 

0.80000E+01 

0.80000E+01 

0.71054E- 14 

0.80000E+01 

0.80000E+01 

0.80000E+01 

0.53291E- 14 

0.80000E+01 

0.80000E+01 

0. 17764E- 14 

0 . 1 7764E - 1 4 

0. 17764E-14 


Table. 8 Element Variations of Li Solution for n = 5 ( 25 Bilinear Elements ) 
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